.libPaths('~/R/x86_64-pc-linux-gnu-library/4.1/')
suppressPackageStartupMessages({
library(plotly)
library(gaston)
library(dplyr)
})
# load engine's functions
sapply(FUN = source,
X = list.files("../../src",
pattern = ".R$",
full.names = T))
# R options
options(max.print = 20)
genoFile <- '../../data/geno/breedGame_phasedGeno.vcf.gz'
crossTableFile <- '../../data/crossingTable/breedGame_small_crossTable.csv'
SNPcoordFile <- '../../data/SNPcoordinates/breedingGame_SNPcoord.csv'
markerEffectsFile <- '../../data/markerEffects/breedGame_markerEffects.csv'
predictions <- calc_progenyBlupEstimation(
genoFile = genoFile,
crossTableFile = crossTableFile,
SNPcoordFile = SNPcoordFile,
markerEffectsFile = markerEffectsFile,
# outFile = outFile
)
predictions$Cross <- paste0(predictions$ind1, '_X_', predictions$ind2)
set.seed(1234)
markerEffects <- readMarkerEffects(markerEffectsFile)
crossSimOutFile <- crossingSimulation(
genoFile = genoFile,
crossTableFile = crossTableFile,
SNPcoordFile = SNPcoordFile,
nCross = 100)
simulation <- gaston::read.vcf(crossSimOutFile)
## Warning in gaston::read.vcf(crossSimOutFile): NAs introduced by coercion
## Warning in gaston::read.vcf(crossSimOutFile): Some unknown chromosomes id's (try
## to set convert.chr = FALSE)
simulation <- gaston::as.matrix(simulation)
simulation <- simulation[, row.names(markerEffects$SNPeffects)]
simulation <- data.frame(
simBV = as.vector(simulation %*% as.matrix(markerEffects$SNPeffects)),
Cross = as.factor(sub('-\\d+$', '', row.names(simulation)))
)
p <- plot_ly()
# add boxplot of simulated BV
p <- p %>% add_boxplot(data = simulation,
x = ~Cross,
y = ~simBV,
name = "Boxplot Simulated BV")
# add jittered markers of simulated BV
p <- p %>% add_boxplot(data = simulation,
x = ~Cross,
y = ~simBV,
line = list(color = 'rgba(0,0,0,0)'),
marker = list(color = 'rgb(31, 119, 180)'),
fillcolor = 'rgba(0,0,0,0)',
name = "Markers Simulated BV",
boxpoints = "all",
pointpos = 0,
jitter = 1,
hoverinfo = 'text',
text = apply(simulation, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
}))
# add predicted BV
p <- p %>% add_markers(
data = predictions,
x = ~Cross,
y = ~blup_exp,
marker = list(color = 'rgb(255,0,0)'),
name = "Predicted BV",
error_y = ~ list(array = sqrt(blup_var),
type = 'data',
color = '#000000'),
hoverinfo = 'text',
text = apply(predictions, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
})
)
p
alpha <- 0.05
beta <- 0.95
pred_vs_sim <- full_join(simulation, predictions, by = "Cross") %>%
filter(Cross != "Coll0659_X_Coll0425") %>%
group_by(Cross) %>%
summarise(obs_mean = mean(simBV),
blup_exp = unique(blup_exp),
p.value = t.test(x = simBV, mu = unique(blup_exp))$p.value,
delta = power.t.test(n = 100, power = beta, sd = sqrt(unique(blup_var)), sig.level = alpha)$delta)
pred_vs_sim$adj.p.val <- p.adjust(pred_vs_sim$p.value, method = "bonferroni")
pred_vs_sim
idLine <- data.frame(mean = c(min(c(pred_vs_sim$blup_exp, pred_vs_sim$obs_mean)) - 1,
max(c(pred_vs_sim$blup_exp, pred_vs_sim$obs_mean)) + 1),
var = c(min(c(pred_vs_sim$blup_var, pred_vs_sim$sim_var)) - 1,
max(c(pred_vs_sim$blup_var, pred_vs_sim$sim_var)) + 1))
## Warning: Unknown or uninitialised column: `blup_var`.
## Warning: Unknown or uninitialised column: `sim_var`.
## Warning in min(c(pred_vs_sim$blup_var, pred_vs_sim$sim_var)): no non-missing
## arguments to min; returning Inf
## Warning: Unknown or uninitialised column: `blup_var`.
## Warning: Unknown or uninitialised column: `sim_var`.
## Warning in max(c(pred_vs_sim$blup_var, pred_vs_sim$sim_var)): no non-missing
## arguments to max; returning -Inf
rmse <- sqrt(mean((pred_vs_sim$blup_exp - pred_vs_sim$obs_mean)^2))
#plot mean
plot_ly(type = "scatter",
mode = "markers",
data = pred_vs_sim,
x = ~blup_exp,
y = ~obs_mean,
name = "Observed mean vs prediction",
hoverinfo = 'text',
text = apply(pred_vs_sim, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
})
) %>%
add_lines(inherit = FALSE,
x = idLine$mean,
y = idLine$mean,
name = "Identity line")
pred_vs_sim <- full_join(simulation, predictions, by = "Cross") %>%
filter(Cross != "Coll0659_X_Coll0425") %>%
group_by(Cross) %>%
summarise(obs_var = var(simBV),
blup_var = unique(blup_var),
p.value = ks.test(simBV, "pnorm", mean = unique(blup_exp), sd=sqrt(unique(blup_var)))$p.value)
pred_vs_sim$adj.p.val <- p.adjust(pred_vs_sim$p.value, method = "bonferroni")
pred_vs_sim
## Document generated in:
## Time difference of 26.87301 secs
##
## CPU: AMD Ryzen 5 3600X 6-Core Processor
## Memory total size: 32.79254 GB
##
##
## Session information:
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Pop!_OS 22.04 LTS
##
## Matrix products: default
## BLAS: /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.0.7 gaston_1.5.7 RcppParallel_5.1.5 Rcpp_1.0.8
## [5] plotly_4.10.0 ggplot2_3.3.5
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.1 xfun_0.29 bslib_0.3.1
## [4] memuse_4.2-1 purrr_0.3.4 splines_4.1.1
## [7] lattice_0.20-44 colorspace_2.0-2 vctrs_0.3.8
## [10] generics_0.1.2 htmltools_0.5.2 viridisLite_0.4.0
## [13] vcfR_1.12.0 mgcv_1.8-36 yaml_2.2.2
## [16] utf8_1.2.2 rlang_1.0.1 jquerylib_0.1.4
## [19] pillar_1.7.0 glue_1.6.1 withr_2.4.3
## [22] DBI_1.1.2 lifecycle_1.0.1 stringr_1.4.0
## [25] munsell_0.5.0 gtable_0.3.0 htmlwidgets_1.5.4
## [28] breedSimulatR_0.3.2 evaluate_0.14 knitr_1.37
## [31] permute_0.9-7 fastmap_1.1.0 crosstalk_1.2.0
## [34] parallel_4.1.1 fansi_1.0.2 pinfsc50_1.2.0
## [37] scales_1.1.1 vegan_2.6-2 jsonlite_1.7.3
## [40] digest_0.6.29 stringi_1.7.6 grid_4.1.1
## [43] cli_3.1.1 tools_4.1.1 magrittr_2.0.2
## [46] sass_0.4.0 lazyeval_0.2.2 tibble_3.1.6
## [49] cluster_2.1.2 ape_5.6-2 crayon_1.4.2
## [52] tidyr_1.2.0 pkgconfig_2.0.3 Matrix_1.3-4
## [55] MASS_7.3-55 ellipsis_0.3.2 data.table_1.14.2
## [58] assertthat_0.2.1 rmarkdown_2.11 httr_1.4.2
## [61] rstudioapi_0.13 R6_2.5.1 nlme_3.1-152
## [64] compiler_4.1.1
shiny::HTML('
<!-- Add icon/font library -->
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.13.0/css/all.min.css">
<link rel="stylesheet" href="https://fonts.googleapis.com/css?family=Lato">
<div class="footer" style="font-family: Lato">
<hr />
<p style="text-align: center;"><a href="https://github.com/juliendiot42">Julien Diot</a></p>
<p style="text-align: center;"><span style="color: #808080;"><em>juliendiot@ut-biomet.org</em></span></p>
<!-- Add font awesome icons -->
<p style="text-align: center;">
<a href="https://github.com/juliendiot42" class="fab fa-github"></a>
<a href="https://www.linkedin.com/in/julien-diot-949592107/" class="fab fa-linkedin"></a>
<a href="https://orcid.org/0000-0002-8738-2034" class="fab fa-orcid"></a>
<a href="https://keybase.io/juliendiot" class="fab fa-keybase"></a>
</p>
</div>
')